
#include "glodef.h"
#include "errlog.h"
#include "grid3d.h"
#include <stdlib.h>
#include <stdio.h>          /* file io */
#include <string.h>        /* memcpy */
// #include <time.h>

#define NEED_GRID_FACES     1
#define NEED_GRID_EDGES     1
#define NEED_NODES_SORTED   1

/****************************************************************************
 *            Generate Grid                                                 *
 ****************************************************************************/

static inline void gerr_alloc(void)
{
    error_msg("GRID3D: Memory allocation error.\n");
}

static void AllocateNodesAndElems(grid, X0, X1, NX, Y0, Y1, NY, Z0, Z1, NZ)
grid3d_t *grid;
double X0, X1, Y0, Y1, Z0, Z1;
int NX, NY, NZ;
{
    node_t *Nd;
    elem_t *El;
    double dX, dY, dZ, X, Y, Z;
    int i, j, k, q, nds[8];
    /* unsigned bc_z, bc_y, bc_x; */
    unsigned zb, yb, xb;

    dX=(X1-X0)/(NX);
    dY=(Y1-Y0)/(NY);
    dZ=(Z1-Z0)/(NZ);

    NX++;
    NY++;
    NZ++;

    grid->NN=NX*NY * NZ;
    grid->NNI=(NX-2)*(NY-2)*(NZ-2);
    /* do nodes first */
    PrintLog("grid3d: Allocating %i nodes... ", grid->NN);
    Nd=(node_t*) calloc(grid->NN, sizeof(node_t));
    if(Nd==NULL) gerr_alloc();
    grid->nodes=Nd;
    for(Z=Z0, k=0; k<NZ; k++, Z+=dZ)
    {
        Nd->on_bdry=0;
        if(k==0) zb=3;
        else if(k==NZ-1) zb=4;
        else zb=0;
        /* bc_z = (k*(k+1-NZ)) ? BC_NON : BC_DIR; */
        for(Y=Y0, j=0; j<NY; j++, Y+=dY)
        {
            if(j==0) yb=2;
            else if(j==NY-1) yb=5;
            else yb=0;
            /* bc_y = (j*(j+1-NZ)) ? bc_z : BC_DIR; */
            for(X=X0, i=0; i<NX; i++, X+=dX)
            {
                if(i==0) xb=1;
                else if(i==NX-1) xb=6;
                else xb=0;
                /* bc_x = (i*(i+1-NZ)) ? bc_y : BC_DIR; */
                Nd->coords[0]=X;
                Nd->coords[1]=Y;
                Nd->coords[2]=Z;
                Nd->on_bdry=(xb?1<<(xb-1):0)|
                        (yb?1<<(yb-1):0)|
                        (zb?1<<(zb-1):0);
                Nd->bcond=(Nd->on_bdry?BC_DIR:BC_NON);
                Nd++;
            }
        }
    }
    PrintLog("done.    ");

    grid->NEL=5*(NX-1)*(NY-1)*(NZ-1);
    /* then do elements */
    PrintLog("Allocating %i elements... ", grid->NEL);
    El=(elem_t*) calloc(grid->NEL, sizeof(elem_t));
    if(El==NULL) gerr_alloc();
    grid->elems=El;
    nds[0]=0;
    nds[1]=nds[0]+1;
    nds[2]=nds[0]+NX;
    nds[3]=nds[2]+1;
    nds[4]=nds[0]+NX*NY;
    nds[5]=nds[4]+1;
    nds[6]=nds[4]+NX;
    nds[7]=nds[6]+1;
    for(k=0; k<NZ-1; k++)
    {
        for(j=0; j<NY-1; j++)
        {
            for(i=0; i<NX-1; i++)
            {
                assert(nds[7]<grid->NN);
                if(nds[0]%2)
                {
                    /* do odd cube */
                    El->nodes[0]=nds[1];
                    El->nodes[1]=nds[2];
                    El->nodes[2]=nds[4];
                    El->nodes[3]=nds[7];
                    El=El+1;
                    El->nodes[0]=nds[0];
                    El->nodes[1]=nds[1];
                    El->nodes[2]=nds[2];
                    El->nodes[3]=nds[4];
                    El=El+1;
                    El->nodes[0]=nds[3];
                    El->nodes[1]=nds[2];
                    El->nodes[2]=nds[1];
                    El->nodes[3]=nds[7];
                    El=El+1;
                    El->nodes[0]=nds[5];
                    El->nodes[1]=nds[1];
                    El->nodes[2]=nds[4];
                    El->nodes[3]=nds[7];
                    El=El+1;
                    El->nodes[0]=nds[6];
                    El->nodes[1]=nds[4];
                    El->nodes[2]=nds[2];
                    El->nodes[3]=nds[7];
                    El=El+1;
                }
                else
                {
                    /* do even cube */
                    El->nodes[0]=nds[5];
                    El->nodes[1]=nds[0];
                    El->nodes[2]=nds[6];
                    El->nodes[3]=nds[3];
                    El=El+1;
                    El->nodes[0]=nds[4];
                    El->nodes[1]=nds[0];
                    El->nodes[2]=nds[6];
                    El->nodes[3]=nds[5];
                    El=El+1;
                    El->nodes[0]=nds[1];
                    El->nodes[1]=nds[3];
                    El->nodes[2]=nds[0];
                    El->nodes[3]=nds[5];
                    El=El+1;
                    El->nodes[0]=nds[2];
                    El->nodes[1]=nds[0];
                    El->nodes[2]=nds[3];
                    El->nodes[3]=nds[6];
                    El=El+1;
                    El->nodes[0]=nds[7];
                    El->nodes[1]=nds[6];
                    El->nodes[2]=nds[3];
                    El->nodes[3]=nds[5];
                    El=El+1;
                }
                assert(El<=grid->elems+grid->NEL);
                for(q=0; q<8; q++) nds[q]+=1;
            }
            for(q=0; q<8; q++) nds[q]+=1;
        }
        for(q=0; q<8; q++) nds[q]+=NX;
    }
    PrintLog("done.\n");

    return;
}

grid3d_t *GenerateGrid(int n)
{
    return GenerateGrid1(0, 1, n, 0, 1, n, 0, 1, n);
}

grid3d_t *GenerateGrid1(X0, X1, NX, Y0, Y1, NY, Z0, Z1, NZ)
double X0, X1, Y0, Y1, Z0, Z1;
int NX, NY, NZ;
{
    grid3d_t *grid;
    char s[1024];

    LOG_MPI_WTIME_START;
    PrintLog("grid3d: GenerateGrid [%g,%g]x[%g,%g]x[%g,%g]@(%ux%ux%u)\n",
             X0, X1, Y0, Y1, Z0, Z1, NX, NY, NZ);
    if((NX%2)==1||(NY%2)==1||(NZ%2)==1)
        error_msg("grid3d: Grid parameter must be even.\n");
    grid=(grid3d_t *) calloc(1, sizeof(grid3d_t));
    if(grid==NULL) gerr_alloc();
    AllocateNodesAndElems(grid, X0, X1, NX, Y0, Y1, NY, Z0, Z1, NZ);
    EnrichGrid(grid);
    UpdateMeasures(grid);
    LogGridStats(grid);
    {
        int n;
        n=sprintf(s, "g%ux%ux%u", NX, NY, NZ)+1;
        grid->gridname=strncpy(calloc((n+1), sizeof(char)), s, n+1);
    }

    LOG_MPI_WTIME_STOP;
    return grid;
}

void LogGridStats(const grid3d_t *grid)
{
#define LOGGIRDSTAT(VAR,DESCR)      \
        if(VAR>0) PrintLog("grid3d: %-45s %i\n",DESCR,VAR)
    LOGGIRDSTAT(grid->NN, "Number of nodes");
    LOGGIRDSTAT(grid->NNI, "Number of nodes inside");
    LOGGIRDSTAT(grid->NE, "Number of edges");
    LOGGIRDSTAT(grid->NEI, "Number of edges inside");
    LOGGIRDSTAT(grid->NDFI, "Number of div-free base functions inside");
    LOGGIRDSTAT(grid->NF, "Number of faces");
    LOGGIRDSTAT(grid->NFI, "Number of faces inside");
    LOGGIRDSTAT(grid->NEL, "Number of elements");
#undef LOGGIRDSTAT
}

void EnrichGrid(grid3d_t *grid)
{
    assert(grid!=NULL);
    assert(grid->nodes!=NULL);
    assert(grid->elems!=NULL);
#ifdef NEED_GRID_FACES
    MakeFaces(grid);
#endif
#ifdef NEED_NODES_SORTED
    SortNodes(grid);
#endif
#ifdef NEED_GRID_EDGES
    MakeEdges(grid);
#endif
}

/****************************************************************************
 *            ReleaseGrid                                                   *
 ****************************************************************************/

void ReleaseGrid(grid3d_t **grid)
{
    PrintLog("grid3d: ReleaseGrid.\n");
    if(*grid==NULL) return;
    if((*grid)->gridname) free((*grid)->gridname);
    if((*grid)->nodes) free((*grid)->nodes);
    if((*grid)->edges) free((*grid)->edges);
    if((*grid)->faces) free((*grid)->faces);
    if((*grid)->elems) free((*grid)->elems);
    if((*grid)->e_faces) free((*grid)->e_faces);

    ReleaseInverse(*grid);

    free(*grid);
    (*grid)=NULL;
}


/****************************************************************************
 *            SortNodes                                                     *
 ****************************************************************************/

static int *OI, *NI;
static node_t *NDS;

static int compNDS_count;

static int compNDS(const void *a, const void *b)
{
    const int *A=a, *B=b;
    int ba, bb;
    compNDS_count++;
    ba=NDS[*A].bcond;
    bb=NDS[*B].bcond;
    return ba-bb;
}

void SortNodes(grid3d_t *grid)
/* OI[i] is the old index of node with new index i */
/* NI[i] is the new index of node with old index i */
{
    int i;
    edge_t *e;
    face_t *f;
    elem_t *T;
    PrintLog("grid3d: SortNodes: ");
    LOG_MPI_WTIME_START;
    /* allocate memory for auxiliary arrays */
    OI=(int*) malloc(grid->NN*sizeof(int));
    NI=(int*) malloc(grid->NN*sizeof(int));
    NDS=(node_t*) malloc(grid->NN*sizeof(node_t));
    if((OI==NULL)||(NI==NULL)||(NDS==NULL)) gerr_alloc();
    /* prepare NDS */
    memcpy(NDS, grid->nodes, grid->NN*sizeof(node_t));
    if(grid->faces)
    {
        for(i=0; i<grid->NN; i++)
            NDS[i].bcond=BC_NON;
        for(f=grid->faces+grid->NFI; f<grid->faces+grid->NF; f++)
            for(i=0; i<3; i++)
                NDS[f->nodes[i]].bcond=BC_DIR;
        grid->NNI=grid->NN;
        for(i=0; i<grid->NN; i++)
            if(NDS[i].bcond==BC_DIR)
                grid->NNI--;
    }
    /* prepare OI -- before sorting OI is identity: */
    for(i=0; i<grid->NN; i++) OI[i]=i;
    /* sort OI */
    compNDS_count=0;
    qsort(OI, grid->NN, sizeof(int), compNDS);
    /* prepare NI -- NI is the inverse of OI: */
    for(i=0; i<grid->NN; i++) NI[OI[i]]=i;
    /* sort grid->nodes */
    for(i=0; i<grid->NN; i++) grid->nodes[i]=NDS[OI[i]];
    /* update node indices in elements, faces, edges */
    if(grid->elems)
        for(T=grid->elems; T<grid->elems+grid->NEL; T++)
            for(i=0; i<4; i++) T->nodes[i]=NI[T->nodes[i]];
    if(grid->faces)
        for(f=grid->faces; f<grid->faces+grid->NF; f++)
            for(i=0; i<3; i++) f->nodes[i]=NI[f->nodes[i]];
    if(grid->edges)
        for(e=grid->edges; e<grid->edges+grid->NE; e++)
            for(i=0; i<2; i++) e->nodes[i]=NI[e->nodes[i]];
    /* done. release memory and leave */
    {
        free(NDS);
        NDS=NULL;
    }
    {
        free(OI);
        OI=NULL;
    }
    {
        free(NI);
        NI=NULL;
    }
    LOG_MPI_WTIME_STOP;
    return;
}

/****************************************************************************
 *            MakeFaces                                                     *
 ****************************************************************************/


typedef struct FL_S
{
    int nodes[3];
    int elems[2];
    int count;
} fl_t;

static int comp_count;

static int comp(const void *a, const void *b)
{
    int A= *((int *) a);
    int B= *((int *) b);
    comp_count++;
    return A-B;
}

static int compFL_count;

static int compFL(const void *A, const void *B)
{
    const int *a=((fl_t*) A)->nodes;
    const int *b=((fl_t*) B)->nodes;
    int c;
    compFL_count++;
    c=(int)*a-(int)*b;
    if(c) return c;
    a++;
    b++;
    c=(int)*a-(int)*b;
    if(c) return c;
    a++;
    b++;
    return(int) *a-(int)*b;
}

static void _set_face(F, n0, n1, n2, el)
fl_t *F;
int n0, n1, n2, el;
{
    int n[3];
    n[0]=n0;
    n[1]=n1;
    n[2]=n2;
    qsort(n, 3, sizeof(int), comp);
    memcpy(F->nodes, n, 3*sizeof(int));
    F->elems[0]=el;
    F->elems[1]= -1;
    F->count=1;
}

void MakeFaces(grid3d_t *grid)
{
    face_t *f;
    fl_t *FF, *F;
    elem_t *T;
    int el, *en, i, j, k, a, b, bc;

    PrintLog("grid3d: MakeFaces: ");
    LOG_MPI_WTIME_START;
    compFL_count=0;
    comp_count=0;
    /* create an array of faces by adding 4 faces for each element */
    F=FF=(fl_t*) malloc(4*grid->NEL*sizeof(fl_t));
    for(el=0, T=grid->elems; el<grid->NEL; el++, T++)
    {
        en=T->nodes;
        _set_face(F++, en[0], en[1], en[2], el);
        _set_face(F++, en[0], en[1], en[3], el);
        _set_face(F++, en[0], en[2], en[3], el);
        _set_face(F++, en[1], en[2], en[3], el);
    }
    /* sort the above array in icreasing order of the indices of the nodes */
    qsort(FF, 4*grid->NEL, sizeof(fl_t), compFL);
    /* combine faces with the same nodes.  inner faces belong to two elements,
       so they will appear twice in the array.  Mark the first occurence with
       count=2 and elems containing both elements.
       Mark the second occurence with count=0 */
    for(F=FF, i=1; i<4*grid->NEL; i++, F++)
    {
        if(compFL(F, FF+i)==0)
        { /* equal */
            F->count=2;
            j=F->elems[0];
            k=FF[i].elems[0];
            minmax(j, k, F->elems, F->elems+1);
            F++;
            F->count=0;
            i++;
        }
    }
    /* count the number of faces and the number of faces inside */
    grid->NF=grid->NFI=0;
    for(F=FF, i=0; i<4*grid->NEL; i++, F++)
    {
        switch(F->count)
        {
            case 2: grid->NFI++; /* do not put "break" here! */
            case 1: grid->NF++;
            case 0: break;
        }
    }
    /* allocate the faces array of the grid */
    grid->faces=(face_t *) malloc(grid->NF*sizeof(face_t));
    if(grid->faces==NULL)
        gerr_alloc();
    /* populate the faces array leaving boundary faces last */
    /* populate element's local .faces arrays as we go */
    a=0;
    b=grid->NFI;
    for(F=FF, i=0; i<4*grid->NEL; i++, F++)
    {
        if(F->count==0) continue;
        if(F->count==2)
        {
            f=grid->faces+a++;
            bc=BC_NON;
        }
        else /* F->count==1 */
        {
            f=grid->faces+b++;
            bc=BC_DIR;
        }
        memcpy(f->nodes, F->nodes, 3*sizeof(int));
        memcpy(f->elems, F->elems, 2*sizeof(int));
        f->bcond=bc;
        for(el=0; el<F->count; el++)
        {
            T=grid->elems+f->elems[el];
            for(j=0; j<4; j++)
            {
                for(k=0; k<3; k++)
                    if(T->nodes[j]==f->nodes[k])
                        break;
                if(k>=3)
                {
                    T->faces[j]=f-grid->faces;
                    break;
                }
            }
            if(j>=4)
                error_msg("MakeFaces: this is impossible.\n");
        }
    }
    free(FF);

    LOG_MPI_WTIME_STOP;
}

/****************************************************************************
 *            MakeEdges                                                     *
 ****************************************************************************/

typedef struct
{
    int nodes[2];
    unsigned bcond;
    int nf;
    int face;
} el_t;

static el_t *EDGS;
static int *NM; /* node marker for creating the tree */

static int compEL_count;

static int compEL(const void *A, const void *B)
{
    const int *a=((el_t*) A)->nodes,
            *b=((el_t*) B)->nodes;
    int c;
    compEL_count++;
    c=(int)*a-(int)*b;
    if(c) return c;
    a++;
    b++;
    return(int) *a-(int)*b;
}

void _set_edge(E, n0, n1, bc, f)
el_t *E;
int n0, n1, bc, f;
{
    minmax(n0, n1, E->nodes, E->nodes+1);
    E->bcond=bc;
    E->nf=1;
    E->face=f;
}

void MakeEdges(grid3d_t *grid)
{
    face_t *f;
    el_t *E;
    elem_t *T;
    int p, m, a, b, c, d;
    edge_t *e;

    if(grid->faces==NULL)
        error_msg("MakeEdges: call MakeFaces first.\n");

    PrintLog("grid3d: MakeEdges: ");
    LOG_MPI_WTIME_START;
    grid->NE=grid->NN+grid->NF-grid->NEL-1;
    /* allocate memory for auxiliary arrays */
    EDGS=(el_t*) calloc(3*grid->NF, sizeof(el_t));
    NM=(int*) calloc(grid->NN, sizeof(int));
    if((EDGS==NULL)||(NM==NULL))
        gerr_alloc();
    /* fill EDGS with 3 edges per face */
    for(E=EDGS, f=grid->faces; f<grid->faces+grid->NF; f++)
    {
        _set_edge(E++, f->nodes[1], f->nodes[2], f->bcond, f-grid->faces);
        _set_edge(E++, f->nodes[0], f->nodes[1], f->bcond, f-grid->faces);
        _set_edge(E++, f->nodes[0], f->nodes[2], f->bcond, f-grid->faces);
    }
    /* sort EDGS in order of increasing node indices */
    compEL_count=0;
    qsort(EDGS, 3*grid->NF, sizeof(el_t), compEL);
    /* after sorting EDGS has groups of edges with the same nodes but different
        ->face.  Mark the first one with ->nf equal to the number of faces on this
        edge and mark the rest with ->nf=0.  Also if any of the faces which added
        this edge is on the boundary then the edge is on the boundary. */
    for(E=EDGS, p=1; p<3*grid->NF; p++)
    {
        if(compEL(E, EDGS+p)==0)
        {
            E->nf++;
            EDGS[p].nf=0;
            if(EDGS[p].bcond!=BC_NON)
                E->bcond=EDGS[p].bcond;
        }
        else
        {
            E=EDGS+p;
        }
    }
    /* Build the spanning tree of the graph of nodes and edges.
        We don't actually store the tree as a data structure, rather, we mark
        the edges in the tree with ->bcond |= BC_ZERO.
        Div-free functions associated with an edge not marked with BC_ZERO
        form a basis of the div-free subspace. */
    E=EDGS;
    NM[grid->NNI]=1;
    m=1;
    while(m<grid->NN)
    {
        int nn, NNB=grid->NN-grid->NNI;
        if(E->nf==0) goto inc_E;
        if(NM[E->nodes[0]]==NM[E->nodes[1]]) goto inc_E;
        if(NM[E->nodes[0]])
            nn=E->nodes[1];
        else
            nn=E->nodes[0];
        if((E->bcond&BC_DIR)||(m>=NNB))
        {
            NM[nn]=1;
            E->bcond|=BC_ZERO;
            m=m+1;
        }
inc_E:
        {
            E++;
            if(E>=EDGS+3*grid->NF)
                E=EDGS;
        }
    }
    /* count the number of edges inside and the number of div-free base
        functions on inside edges.
        We have the following situation:

        0                inside edges            NEI    bdry edges         NE
        |--------------------------------|--------|---------|---------------|
                                       NDFI
                                         |  marked BC_ZERO  |
                                         |      NN - 1      |
     */
    grid->NEI=0;
    grid->NDFI=0;
    for(E=EDGS, p=0; p<3*grid->NF; p++, E++)
    {
        if(E->nf==0)
            continue;
        if((E->bcond&BC_DIR)==0)
        {
            grid->NEI++;
            if((E->bcond&BC_ZERO)==0)
                grid->NDFI++;
        }
    }
    /* allocate memory */
    grid->e_faces=(int*) malloc(3*grid->NF*sizeof(int));
    grid->edges=(edge_t*) malloc(grid->NE*sizeof(edge_t));
    if((grid->edges==NULL)||(grid->e_faces==NULL)) gerr_alloc();
    /* populate the e_faces array */
    for(m=0; m<3*grid->NF; m++) grid->e_faces[m]=EDGS[m].face;
    /* populate the edges array */
    a=0;
    b=grid->NDFI;
    c=grid->NEI;
    d=grid->NDFI+grid->NN-1;
    for(E=EDGS, p=0; p<3*grid->NF; p++, E++)
    {
        if(E->nf==0) continue;
        if((E->bcond&BC_DIR)&&(E->bcond&BC_ZERO))
            e=grid->edges+c++;
        else if(E->bcond&BC_DIR)
            e=grid->edges+d++;
        else if(E->bcond&BC_ZERO)
            e=grid->edges+b++;
        else
            e=grid->edges+a++;
        e->nodes[0]=E->nodes[0];
        e->nodes[1]=E->nodes[1];
        e->bcond=E->bcond;
        e->nf=E->nf;
        e->faces=grid->e_faces+p;
    }
    /* update the edges arrays in faces */
    for(e=grid->edges, p=0; p<grid->NE; p++, e++)
    {
        for(m=0; m<e->nf; m++)
        {
            f=grid->faces+e->faces[m];
            for(a=0; a<3; a++)
            {
                if((f->nodes[a]==e->nodes[0])||(f->nodes[a]==e->nodes[1]))
                    continue;
                /* else */
                f->edges[a]=p;
                break;
            }
        }
    }
    /* update the edges arrays in elems */
    for(T=grid->elems; T<grid->elems+grid->NEL; T++)
    {
        f=grid->faces+T->faces[3];
        for(m=0; m<3; m++)
        {
            if(f->nodes[m]==T->nodes[0])
                T->edges[3]=f->edges[m];
            else
                if(f->nodes[m]==T->nodes[1])
                T->edges[1]=f->edges[m];
            else /* if(f->nodes[m]==T->nodes[2]) */
                T->edges[0]=f->edges[m];
        }

        f=grid->faces+T->faces[1];
        for(m=0; m<3; m++)
        {
            if(f->nodes[m]==T->nodes[0])
                T->edges[5]=f->edges[m];
            else
                if(f->nodes[m]==T->nodes[2])
                T->edges[2]=f->edges[m];
        }

        f=grid->faces+T->faces[0];
        for(m=0; m<3; m++)
        {
            if(f->nodes[m]==T->nodes[2])
            {
                T->edges[4]=f->edges[m];
                break;
            }
        }
    }

    /* clean up and leave */
    free(EDGS);
    free(NM);

    LOG_MPI_WTIME_STOP;
}

/****************************************************************************
 *                Inverse connectvity                                       *
 ****************************************************************************/

void MakeInverse(grid3d_t *grid)
{
    inverse_t *inv;
    node_inv_t *ninv;
    /* edge_inv_t *einv; */
    void *V;
    int id;
    edge_t *e;
    face_t *f;
    elem_t *T;

    PrintLog("grid3d: MakeInverse: ");
    LOG_MPI_WTIME_START;

    if(grid==NULL) return;
    if(grid->inverse)
        ReleaseInverse(grid);
    /* allocate memory for inverse connectivity table */
    V=calloc(1, sizeof(inverse_t));
    if(V==NULL)
        gerr_alloc();
    grid->inverse=inv=V;
    /* step 1: inverse connectivity for nodes */
    /* allocate memory for inverse connectivity for nodes */
    V=calloc(grid->NN, sizeof(node_inv_t));
    if(V==NULL)
        gerr_alloc();
    inv->ninv=ninv=V;
    /* step 1.1: fill in connected nodes first */
    /* count the number of connected nodes (== number of containing edges) */
    for(e=grid->edges, id=0; id<grid->NE; id++, e++)
    {
        ninv[e->nodes[0]].NN++;
        ninv[e->nodes[1]].NN++;
    }
    /* allocate memory for local .nodes arrays. clean up .NN */
    for(id=0; id<grid->NN; id++)
    {
        V=calloc(ninv[id].NN, sizeof(int));
        if(V==NULL)
            gerr_alloc();
        ninv[id].nids=V;
        ninv[id].NN=0;
    }
    /* fill in the local .nids arrays */
    for(e=grid->edges, id=0; id<grid->NE; id++, e++)
    {
        ninv[e->nodes[0]].nids[ninv[e->nodes[0]].NN++]=e->nodes[1];
        ninv[e->nodes[1]].nids[ninv[e->nodes[1]].NN++]=e->nodes[0];
    }
    /* step 1.2: now fill in the containing elements */
    /* count the number of containing elements */
    for(T=grid->elems, id=0; id<grid->NEL; id++, T++)
    {
        ninv[T->nodes[0]].NEL++;
        ninv[T->nodes[1]].NEL++;
        ninv[T->nodes[2]].NEL++;
        ninv[T->nodes[3]].NEL++;
    }
    /* allocate memory for local .elems arrays. clean up .NEL */
    for(id=0; id<grid->NN; id++)
    {
        V=calloc(ninv[id].NEL, sizeof(int));
        if(V==NULL)
            gerr_alloc();
        ninv[id].Tids=V;
        ninv[id].NEL=0;
    }
    /* now populate the .elems arrays */
    for(T=grid->elems, id=0; id<grid->NEL; id++, T++)
    {
        ninv[T->nodes[0]].Tids[ninv[T->nodes[0]].NEL++]=id;
        ninv[T->nodes[1]].Tids[ninv[T->nodes[1]].NEL++]=id;
        ninv[T->nodes[2]].Tids[ninv[T->nodes[2]].NEL++]=id;
        ninv[T->nodes[3]].Tids[ninv[T->nodes[3]].NEL++]=id;
    }
    /* step 1.3: now fill in the containing faces */
    /* count the number of containing faces */
    for(f=grid->faces, id=0; id<grid->NF; id++, f++)
    {
        ninv[f->nodes[0]].NF++;
        ninv[f->nodes[1]].NF++;
        ninv[f->nodes[2]].NF++;
    }
    /* allocate memory for local .fids arrays. clean up .NF */
    for(id=0; id<grid->NN; id++)
    {
        V=calloc(ninv[id].NF, sizeof(int));
        if(V==NULL)
            gerr_alloc();
        ninv[id].fids=V;
        ninv[id].NF=0;
    }
    /* now populate the .fids arrays */
    for(f=grid->faces, id=0; id<grid->NF; id++, f++)
    {
        ninv[f->nodes[0]].fids[ninv[f->nodes[0]].NF++]=id;
        ninv[f->nodes[1]].fids[ninv[f->nodes[1]].NF++]=id;
        ninv[f->nodes[2]].fids[ninv[f->nodes[2]].NF++]=id;
    }
    /* step 2: do inverse connectivity for edges */
    /* allocate memory for inverse connectivity for nodes */
    /*
        V = calloc(grid->NE, sizeof(edge_inv_t));
        if(V==NULL)
            gerr_alloc();
        inv->einv = einv = V;
     */
    /* step 2.1: fill in connected nodes first */
    /* count the number of connected nodes (== number of containing faces */
    /*
        for(f=grid->faces,id=0;id<grid->NF;id++,f++)
        {
            einv[f->edges[0]].NN++;
            einv[f->edges[1]].NN++;
            einv[f->edges[2]].NN++;
        }
     */
    /* allocate memory for .nodes arrays */
    /*
        for(id=0;id<grid->NE;id++)
        {
            V = calloc(einv[id].NN,sizeof(int));
            if(V==NULL)
                gerr_alloc();
            einv[id].nids = V;
            einv[id].NN = 0;
        }
     */
    /* now populate the .elems arrays */
    /*
        for(f=grid->faces,id=0;id<grid->NF;id++,f++)
        {
            einv[f->edges[0]].nids[einv[f->edges[0]].NN++] = f->nodes[0];
            einv[f->edges[1]].nids[einv[f->edges[1]].NN++] = f->nodes[1];
            einv[f->edges[2]].nids[einv[f->edges[2]].NN++] = f->nodes[2];
        }
     */
    /* done :-) */
    LOG_MPI_WTIME_STOP;
    return;
}

void ReleaseInverse(grid3d_t *grid)
{
    if(grid->inverse)
    {
        int id;
        if(grid->inverse->ninv)
        {
            for(id=0; id<grid->NN; id++)
            {
                if(grid->inverse->ninv[id].nids)
                    free(grid->inverse->ninv[id].nids);
                if(grid->inverse->ninv[id].Tids)
                    free(grid->inverse->ninv[id].Tids);
            }
            free(grid->inverse->ninv);
        }
        /*
        if(grid->inverse->einv)
        {
            for(id=0;id<grid->NE;id++)
            {
                if(grid->inverse->einv[id].nids)
                    free(grid->inverse->einv[id].nids);
            }
            free(grid->inverse->einv);
        }
         */
        free(grid->inverse);
        grid->inverse=NULL;
    }
}
